/*
 ============================================================================
 Name        : mBuild_inhomogenous.cpp
 Author      : Yi He
 Version     :
 Copyright   : A free distributable Academic program
 Description : The main program to build a membrane.
 ============================================================================
 */

#include <iostream>
#include "time.h"
#include "mbuild.h"

using namespace std;
void find_file(char* inFile,char* type,char* readfile,char* outfile);
void build_diel(char* fname,ToroidalPore &pore);
void build_charge(char* fname,ToroidalPore &pore);
void build_kappa(char* fname,ToroidalPore &pore);
void build_phi(char* fname,ToroidalPore &pore);
void build_shape(char* fname,ToroidalPore &pore);
void setup_pore(ToroidalPore&);

int main(void) {
	/* Input a file name to guess the other files*/

	  char inFile[256];
	  cout<<"This program build the membrane solvated protein system for Poisson-Boltzmann calculation with APBS.";
	  cout<<"It reads in the DX maps for water solvated protein system and generates DX maps readable for APBS.";
	  cout<<"For DX file standard, please refer to APBS webpage.";
      cout<<"The input files should follow the nameing standard of: \n";
      cout<<"dielx_X.dx,diely_X.dx,dielz_X.dx --- x/y/z-shifted dielectric map.\n";
      cout<<"charge_X.dx --- charge density map.\n";
      cout<<"kappa_X.dx --- ion permibility map.\n";
      cout<<"with the character X in the above files meaning the focusing level.\n";
      cout<<"The above files can be easily generated from dummy run of APBS.";
      cout<<"The result DX maps will have names as dielx_Xm.dx,charge_Xm.dx ...\n";
      cout<<"===================================================================\n";
	  cout<<"Input a file name to get all other files:\n";
	  cin>>inFile;
	  //Setup pore
	  ToroidalPore pore;
	  setup_pore(pore);
	  build_diel(inFile,pore);
	  build_charge(inFile,pore);
	  build_kappa(inFile,pore);
	  build_phi(inFile,pore);
	  build_shape(inFile,pore);
	  return 0;
}

void setup_pore(ToroidalPore &pore)
{
  double r0,k0,xl,t,xh,yl,yh,zl,zh,cutoff;
  int ncyc;
  double r_top,r_bottom;
  cout<<"Set up toroidal pore size(set R0=0 and K0=0 if don't want to have a toroidal pore).";
  cout<<"In current model, torodial pore is considered to be a curved surface following equation:\n\n";
  cout<<"r= R0+ K0*(2z/T)^2,\n\n in which r,z are the radial and axial coordinates, T is the membrane thicknes.";
  cout<<" There are many definition of membrane thicknes but in our case it is the hydrocarbon core thickness.\n";
  cout<<"===================================================================\n";
  cout<<"Pore radius R0 = ";
  cin>>r0;
  cout<<"Curvature factor K0 =";
  cin>>k0;
  cout<<"Membrane thickness T=";
  cin>>t;
  cout<<"Set up hole(for cylindrical pore)\n";
  cout<<"The incontinuity of cylindrical pore is created by creating a hole filled with water inside hydrophobic core of membrane.";
  cout<<"In a cylindrical(or cone shaped) region which have R_bottom at the bottom leaflet and R_top at the top leaflet, the properties(dielectric constant, charge density...etc.) are swithed to those of water.\n";
  cout<<"===================================================================\n";
  cout<<"R_bottom=";
  cin>>r_bottom;
  cout<<"R_top=";
  cin>>r_top;
  cout<<"Build region:\n";
  cout<<"Inside a defined box, the values of input box will be replaced with that of a membrane system.";
  cout<<"If the defined box is larger than the input box, the whole input box will be filled with membrane system.\n";
  cout<<"Please input the range of X,Y,Z for the membrane constructing box in the following:\n";
  cout<<"===================================================================\n";
  cout<<"X min=";
  cin>>xl;
  cout<<"X max=";
  cin>>xh;
  cout<<"Y min=";
  cin>>yl;
  cout<<"Y max=";
  cin>>yh;
  cout<<"Z min=";
  cin>>zl;
  cout<<"Z max=";
  cin>>zh;
  cout<<"Setup integration parameter:\n";
  cout<<"For accurately reproduce the values in each grid points, volume integration is used for points close to membrane. ";
  cout<<"The cutoff defines the range of grid points that will use integration over the grid volume instead of single point value in the grid center. ";
  cout<<"The integration steps defines the how many integration iteration will be used. ";
  cout<<"The longer integration cutoff and the larger integration steps, the more accurate the values will although longer time will be used. ";
  cout<<"10 and 3-6 could be good enough for preparing systems for PB calculations.\n";
  cout<<"===================================================================\n";
  cout<<"Integration cutoff:";
  cin>>cutoff;
  cout<<"Integration steps:";
  cin>>ncyc;
  pore.set_region(xl,xh,yl,yh,zl,zh);
  pore.set_cutoff(cutoff);
  pore.set_integration_cyc(ncyc);
  pore.setup(r0,k0,t);
  CylindricalHole hole=CylindricalHole(r_top,r_bottom,t+10*cutoff,-t/2-5*cutoff);
  if (r_top >0 && r_bottom >0) {
      pore.has_hole(&hole);
   }
}

void build_diel(char* fname,ToroidalPore &pore)
{
  double memb,head,water,width,diel_p;
  //Setup Switch for dielectric constant;
  cout<<"Set up the dielectric distribution:\n";
  cout<<"In the input DX maps, the solute(protein or nucleic acid) may be constructed with a dielectric constant. ";
  cout<<"You should specify that value so that the solute occupied grids will be recognized thus not be overrided.\n";
  cout<<"===================================================================\n";
  cout<<"Dielectric constant in the hydrocarbon core of the membrane:";
  cin>>memb;
  cout<<"Dielectric constant in Head group region:";
  cin>>head;
  cout<<"Dielectric constant in Water:";
  cin>>water;
  cout<<"Head group region  width:";
  cin>>width;
  cout<<"Solute dielectric constant(to skip):";
  cin>>diel_p;
  DielectricConstant c;
  c.set_diels(water,head,memb,diel_p);
  c.set_headgroup_width(width);
  pore.set_distribution(&c);
  GridSpace grid,reference;
  char input[256];
  char output[256];
  find_file(fname,"dielx",input,output);
  cout<<"Reading "<<input<<endl;
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);
  find_file(fname,"diely",input,output);
  cout<<"Reading "<<input<<endl;
  t1=clock();
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
  t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  find_file(fname,"dielz",input,output);
  cout<<"Reading "<<input<<endl;
  t1=clock();
  DxFile::read_dx(input,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Building "<<output<<" ...\n";
  pore.fill(&grid,&reference);
  cout<<"Creating hole\n";
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
  t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
}

void build_charge(char *fname,ToroidalPore &pore)
{
  double offset,area,anfr,width;
  cout<<"Setup charge distribution.\n";
  GridSpace grid,reference;
  int n=2;
  char read[256];
  char output[256];
  ChargeDensity c=ChargeDensity(n);
  QuadraticChargeDistribution d;
  double a,b;
  find_file(fname,"charge",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  cout<<"Head group is Homogeneous on the toroidal pore surface? If not, what is the percentage of charge in the center compared to the pore rim?";
  cout<<"In our model, the charge density is considered to change quadratically with z if the membrane is not homogeneous. Notice that this is only effective for toroidal pore.\n";
  cout<<"Input 1 (lipid density at pore center is 100% that of planar membrane) if it is homogeneous, 0.0 to ? for the type of inhomogeneous pore you have.\n";
  cout<<"===================================================================\n";
  cout<<"Input(0.0-1.0):";
  cin>>a;
  //a=(20*pore.r0*pore.k0/pore.T+6*pore.k0*pore.k0/pore.T-5*pore.r0-3*pore.k0)/(10*pore.r0+2*pore.k0);
  b=4*(1-a)/(pore.T*pore.T);
  d.setup(a,b,pore.r0,pore.k0,pore.T);
  cout<<"Head group has N layers of charge? The more charge layers you used to mimic the charge distribution in head group region, the more close to realistic your model will be, even though computing time will increase dramatically with increasing N.\n";
  cout<<"===================================================================\n";
  cout<<"N = ";
  cin>>n;
  c.set_nlayer(n);
  cout<<"Input charge density for each charge layer:\n";
  cout<<"Has "<<n<<" layers.\n";
  for(int i=1; i<=n;i++)
  {
  cout<<"Layer  "<< i<<":\n";
  cout<<"===================================================================\n";
  cout<<"Charge per lipid:";
  cin>>anfr;
  cout<<"Area per lipid:";
  cin>>area;
  cout<<"Offset of charged plane (from surface of hydrocarbon core):";
  cin>>offset;
  cout<<"Gaussian distribution width:";
  cin>>width;
  c.set_layer(i-1,offset,width,anfr/area);
  }
  cout<<"Building "<<output<<" ...\n";
  pore.set_distribution(&c,&d);
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  pore.create_hole(&grid,&reference);
  DxFile::write_dx(output,&grid);
}

void build_kappa(char* fname,ToroidalPore &pore)
{
  double width;
  cout<<"Setup kappa.\n";
  cout<<"Set ion inacessible region width (from hydrocarbon core, hydrocarbon core is always inpenetrable for ions):\n";
  cout<<"===================================================================\n";
  cout<<"width = ";
  cin>>width;
  Kappa k;
  k.set_headgroup_width(width);
  pore.set_distribution(&k);
  GridSpace grid,reference;
  char read[256];
  char output[256];
  find_file(fname,"kappa",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  grid.copyto(&reference);
  cout<<"Finished.\n";
  clock_t t1=clock();
  pore.fill(&grid,&reference);
  pore.create_hole(&grid,&reference);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);

}

void build_phi(char* inFile,ToroidalPore &pore)
{
  Gouy g;
  double temp,offset,conc,z,area,anfr;
  cout<<"Generating electric potential based on Gouy-Chapmann theory.\n";
  cout<<"===================================================================\n";
  cout<<"Temperature =";
  cin>>temp;
  g.set_temp(temp);
  cout<<"Charge plane offset (from hydrocarbon core) =";
  cin>>offset;
  g.set_offset(offset);
  cout<<"Ion concentration =";
  cin>>conc;
  cout<<"Valence of ion =";
  cin>>z;
  g.set_ion(z,conc);
  cout<<"Area per Lipid =";
  cin>>area;
  cout<<"Fraction of anionic lipids =";
  cin>>anfr;
  g.set_lipid(area,anfr,1.0);
  pore.set_distribution(&g);
  GridSpace grid;
  char read[256];
  char output[256];
  find_file(inFile,"charge",read,output);
  cout<<"Reading "<<read<<"..."<<endl;
  DxFile::read_dx(read,&grid);
  cout<<"Finished.\n";
  grid.clear(0);
  find_file(inFile,"phi",read,output);
  cout<<"Building "<<output<<" ...\n";
  clock_t t1=clock();
  pore.fill(&grid,&grid);
  pore.create_hole(&grid,&grid);
  clock_t t2=clock();
  cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
  DxFile::write_dx(output,&grid);

}
void build_shape(char* fname,ToroidalPore &pore)
{
   GridSpace grid;
   char read[256];
   char output[256];
   IMM1Transition shape;
   double R0,K0,T;
   R0=pore.r0;
   K0=pore.k0;
   T=pore.T;
   //DEBUG
//   cout<<"Pore "<<R0<<" "<<K0<<" "<<T;
   shape.set_nsmth(10);
   shape.set_memb_thickness(T);
   pore.set_distribution(&shape);
   find_file(fname,"charge",read,output);
   DxFile::read_dx(read,&grid);
   grid.clear(0);
   find_file(fname,"shape",read,output);
   cout<<"Building "<<read <<endl;
   clock_t t1=clock();
   pore.fill(&grid,&grid);
   pore.create_hole(&grid,&grid);
   clock_t t2=clock();
   cout<<"Built in "<<double(t2-t1)/CLOCKS_PER_SEC<<"seconds.\n";
   DxFile::write_dx(output,&grid);

}
void find_file(char* inFile,char* type,char* readfile,char* outfile)
{
  int l=strlen(inFile);
  char mid[256];
  strncpy(mid,&inFile[5],l-8);
  int l1=strlen(type);
  strncpy(readfile,type,l1);
  readfile[l1]='\0';
  strcat(readfile,mid);
  readfile[l+l1-8]='\0';
  strcpy(outfile,readfile);
  strcat(readfile,".dx");
  outfile[l+l1-8]='\0';
  strcat(outfile,"m.dx");
}
